<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>dsterf.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="DSTERF.1"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a>( N, D, E, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      INTEGER            INFO, N
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      DOUBLE PRECISION   D( * ), E( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="DSTERF.17"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a> computes all eigenvalues of a symmetric tridiagonal matrix
</span><span class="comment">*</span><span class="comment">  using the Pal-Walker-Kahan variant of the QL or QR algorithm.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (input/output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, the n diagonal elements of the tridiagonal matrix.
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, the eigenvalues in ascending order.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          On entry, the (n-1) subdiagonal elements of the tridiagonal
</span><span class="comment">*</span><span class="comment">          matrix.
</span><span class="comment">*</span><span class="comment">          On exit, E has been destroyed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          &gt; 0:  the algorithm failed to find all of the eigenvalues in
</span><span class="comment">*</span><span class="comment">                a total of 30*N iterations; if INFO = i, then i
</span><span class="comment">*</span><span class="comment">                elements of E have not converged to zero.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE, TWO, THREE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
     $                   THREE = 3.0D0 )
      INTEGER            MAXIT
      PARAMETER          ( MAXIT = 30 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
     $                   NMAXIT
      DOUBLE PRECISION   ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
     $                   OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
     $                   SIGMA, SSFMAX, SSFMIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMCH.59"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.59"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>, <a name="DLAPY2.59"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>
      EXTERNAL           <a name="DLAMCH.60"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.60"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>, <a name="DLAPY2.60"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DLAE2.63"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a>, <a name="DLASCL.63"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>, <a name="DLASRT.63"></a><a href="dlasrt.f.html#DLASRT.1">DLASRT</a>, <a name="XERBLA.63"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, SIGN, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.LT.0 ) THEN
         INFO = -1
         CALL <a name="XERBLA.78"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DSTERF.78"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a>'</span>, -INFO )
         RETURN
      END IF
      IF( N.LE.1 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine the unit roundoff for this environment.
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.86"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'E'</span> )
      EPS2 = EPS**2
      SAFMIN = <a name="DLAMCH.88"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> )
      SAFMAX = ONE / SAFMIN
      SSFMAX = SQRT( SAFMAX ) / THREE
      SSFMIN = SQRT( SAFMIN ) / EPS2
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the eigenvalues of the tridiagonal matrix.
</span><span class="comment">*</span><span class="comment">
</span>      NMAXIT = N*MAXIT
      SIGMA = ZERO
      JTOT = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine where the matrix splits and choose QL or QR iteration
</span><span class="comment">*</span><span class="comment">     for each block, according to whether top or bottom diagonal
</span><span class="comment">*</span><span class="comment">     element is smaller.
</span><span class="comment">*</span><span class="comment">
</span>      L1 = 1
<span class="comment">*</span><span class="comment">
</span>   10 CONTINUE
      IF( L1.GT.N )
     $   GO TO 170
      IF( L1.GT.1 )
     $   E( L1-1 ) = ZERO
      DO 20 M = L1, N - 1
         IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
     $       1 ) ) ) )*EPS ) THEN
            E( M ) = ZERO
            GO TO 30
         END IF
   20 CONTINUE
      M = N
<span class="comment">*</span><span class="comment">
</span>   30 CONTINUE
      L = L1
      LSV = L
      LEND = M
      LENDSV = LEND
      L1 = M + 1
      IF( LEND.EQ.L )
     $   GO TO 10
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale submatrix in rows and columns L to LEND
</span><span class="comment">*</span><span class="comment">
</span>      ANORM = <a name="DLANST.130"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>( <span class="string">'I'</span>, LEND-L+1, D( L ), E( L ) )
      ISCALE = 0
      IF( ANORM.GT.SSFMAX ) THEN
         ISCALE = 1
         CALL <a name="DLASCL.134"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
     $                INFO )
         CALL <a name="DLASCL.136"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
     $                INFO )
      ELSE IF( ANORM.LT.SSFMIN ) THEN
         ISCALE = 2
         CALL <a name="DLASCL.140"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
     $                INFO )
         CALL <a name="DLASCL.142"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
     $                INFO )
      END IF
<span class="comment">*</span><span class="comment">
</span>      DO 40 I = L, LEND - 1
         E( I ) = E( I )**2
   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Choose between QL and QR iteration
</span><span class="comment">*</span><span class="comment">
</span>      IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
         LEND = LSV
         L = LENDSV
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( LEND.GE.L ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        QL Iteration
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Look for small subdiagonal element.
</span><span class="comment">*</span><span class="comment">
</span>   50    CONTINUE
         IF( L.NE.LEND ) THEN
            DO 60 M = L, LEND - 1
               IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
     $            GO TO 70
   60       CONTINUE
         END IF
         M = LEND
<span class="comment">*</span><span class="comment">
</span>   70    CONTINUE
         IF( M.LT.LEND )
     $      E( M ) = ZERO
         P = D( L )
         IF( M.EQ.L )
     $      GO TO 90
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If remaining matrix is 2 by 2, use <a name="DLAE2.179"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a> to compute its
</span><span class="comment">*</span><span class="comment">        eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.EQ.L+1 ) THEN
            RTE = SQRT( E( L ) )
            CALL <a name="DLAE2.184"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a>( D( L ), RTE, D( L+1 ), RT1, RT2 )
            D( L ) = RT1
            D( L+1 ) = RT2
            E( L ) = ZERO
            L = L + 2
            IF( L.LE.LEND )
     $         GO TO 50
            GO TO 150
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( JTOT.EQ.NMAXIT )
     $      GO TO 150
         JTOT = JTOT + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Form shift.
</span><span class="comment">*</span><span class="comment">
</span>         RTE = SQRT( E( L ) )
         SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
         R = <a name="DLAPY2.202"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( SIGMA, ONE )
         SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
<span class="comment">*</span><span class="comment">
</span>         C = ONE
         S = ZERO
         GAMMA = D( M ) - SIGMA
         P = GAMMA*GAMMA
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Inner loop
</span><span class="comment">*</span><span class="comment">
</span>         DO 80 I = M - 1, L, -1
            BB = E( I )
            R = P + BB
            IF( I.NE.M-1 )
     $         E( I+1 ) = S*R
            OLDC = C
            C = P / R
            S = BB / R
            OLDGAM = GAMMA
            ALPHA = D( I )
            GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
            D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
            IF( C.NE.ZERO ) THEN
               P = ( GAMMA*GAMMA ) / C
            ELSE
               P = OLDC*BB
            END IF
   80    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         E( L ) = S*P
         D( L ) = SIGMA + GAMMA
         GO TO 50
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Eigenvalue found.
</span><span class="comment">*</span><span class="comment">
</span>   90    CONTINUE
         D( L ) = P
<span class="comment">*</span><span class="comment">
</span>         L = L + 1
         IF( L.LE.LEND )
     $      GO TO 50
         GO TO 150
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        QR Iteration
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Look for small superdiagonal element.
</span><span class="comment">*</span><span class="comment">
</span>  100    CONTINUE
         DO 110 M = L, LEND + 1, -1
            IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
     $         GO TO 120
  110    CONTINUE
         M = LEND
<span class="comment">*</span><span class="comment">
</span>  120    CONTINUE
         IF( M.GT.LEND )
     $      E( M-1 ) = ZERO
         P = D( L )
         IF( M.EQ.L )
     $      GO TO 140
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If remaining matrix is 2 by 2, use <a name="DLAE2.265"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a> to compute its
</span><span class="comment">*</span><span class="comment">        eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.EQ.L-1 ) THEN
            RTE = SQRT( E( L-1 ) )
            CALL <a name="DLAE2.270"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a>( D( L ), RTE, D( L-1 ), RT1, RT2 )
            D( L ) = RT1
            D( L-1 ) = RT2
            E( L-1 ) = ZERO
            L = L - 2
            IF( L.GE.LEND )
     $         GO TO 100
            GO TO 150
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( JTOT.EQ.NMAXIT )
     $      GO TO 150
         JTOT = JTOT + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Form shift.
</span><span class="comment">*</span><span class="comment">
</span>         RTE = SQRT( E( L-1 ) )
         SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
         R = <a name="DLAPY2.288"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( SIGMA, ONE )
         SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
<span class="comment">*</span><span class="comment">
</span>         C = ONE
         S = ZERO
         GAMMA = D( M ) - SIGMA
         P = GAMMA*GAMMA
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Inner loop
</span><span class="comment">*</span><span class="comment">
</span>         DO 130 I = M, L - 1
            BB = E( I )
            R = P + BB
            IF( I.NE.M )
     $         E( I-1 ) = S*R
            OLDC = C
            C = P / R
            S = BB / R
            OLDGAM = GAMMA
            ALPHA = D( I+1 )
            GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
            D( I ) = OLDGAM + ( ALPHA-GAMMA )
            IF( C.NE.ZERO ) THEN
               P = ( GAMMA*GAMMA ) / C
            ELSE
               P = OLDC*BB
            END IF
  130    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         E( L-1 ) = S*P
         D( L ) = SIGMA + GAMMA
         GO TO 100
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Eigenvalue found.
</span><span class="comment">*</span><span class="comment">
</span>  140    CONTINUE
         D( L ) = P
<span class="comment">*</span><span class="comment">
</span>         L = L - 1
         IF( L.GE.LEND )
     $      GO TO 100
         GO TO 150
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Undo scaling if necessary
</span><span class="comment">*</span><span class="comment">
</span>  150 CONTINUE
      IF( ISCALE.EQ.1 )
     $   CALL <a name="DLASCL.337"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
     $                D( LSV ), N, INFO )
      IF( ISCALE.EQ.2 )
     $   CALL <a name="DLASCL.340"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
     $                D( LSV ), N, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check for no convergence to an eigenvalue after a total
</span><span class="comment">*</span><span class="comment">     of N*MAXIT iterations.
</span><span class="comment">*</span><span class="comment">
</span>      IF( JTOT.LT.NMAXIT )
     $   GO TO 10
      DO 160 I = 1, N - 1
         IF( E( I ).NE.ZERO )
     $      INFO = INFO + 1
  160 CONTINUE
      GO TO 180
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Sort eigenvalues in increasing order.
</span><span class="comment">*</span><span class="comment">
</span>  170 CONTINUE
      CALL <a name="DLASRT.357"></a><a href="dlasrt.f.html#DLASRT.1">DLASRT</a>( <span class="string">'I'</span>, N, D, INFO )
<span class="comment">*</span><span class="comment">
</span>  180 CONTINUE
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DSTERF.362"></a><a href="dsterf.f.html#DSTERF.1">DSTERF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
